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ABSTRACT 

In  order  to  study  multj-dimensional  unstable  detonation  waves,  we  have  developed  a  high 
order  numerical  scheme  suitable  for  calculating  the  detailed  transverse  wave  structures  of 
multidimensional  detonation  waves.  The  numerical  algorithm  uses  a  multi-domain  approach 
so  different  numerical  techniques  can  be  applied  for  different  components  of  detonation  waves. 
The  detonation  waves  are  assumed  to  undergo  an  irreversible,  unimolecular  reaction  A  — ► 
B.  Several  cases  of  unstable  two  dimensional  detonation  waves  are  simulated  and  detailed 
transverse  wave  interactions  are  documented.  The  numerical  results  show  the  importance  of 
resolving  the  detonation  front  without  excessive  numerical  viscosity  in  order  to  obtain  the 
correct  cellular  patterns. 
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Introduction 


Detonation  waves  are  intrinsically  multi-dimensional  unstable  phenomenon  as  demon¬ 
strated  evidently  by  the  experiments  of  Oppenheim  [1]  and  White  [2].  Since  then,  the  simple 
steady  Chapman-Jouquet  theory  [3],  [4]  has  been  re-examined  for  its  limited  explaination 
of  most  of  the  multi-dimensional  features  seen  in  real  world  detonation  waves.  The  unique 
characteristic  of  multi-dimensional  detonation  waves  are  their  cellular  patterns  which  are 
the  trajectories  recorded  on  the  wall  by  the  transverse  waves  structures.  Those  transverse 
wave  structures  consist  of  so-called  “triple  points”  which  are  of  three  shock  configurations 
(an  incident  shock,  a  reflected  shock  and  a  Mach  stem  plus  a  contact  discontinuity)  [5].  The 
formation  of  those  transverse  wave  structures  moving  along  the  main  precursor  detonation 
front  attracted  much  attention.  These  phenomena  raised  the  interest  of  experimentalists 
trying  to  measure  the  cell  size  of  the  cellular  pattern  [6]. 

There  are  many  aspects  in  the  investigations  of  detonation  waves.  Most  important  is 
the  formation  of  a  compressible  detonation  wave  from  laminar  deflagration  waves  [7], [8] 
-  a  process  called  the  DDT  (Deflagration  to  Detonation  Transition)  process.  In  a  DDT 
process,  turbulent  boundary  layers  are  recognized  as  playing  an  important  role  in  flame 
acceleration  and  the  buildup  of  compressible  pressure  waves  in  front  of  flame  fronts.  The 
latter  eventually  will  cause  ‘explosions  in  explosions’  in  the  reacting  flows  and  generate 
detonation  waves.  Applied  mathematicians  have  made  different  attempts  to  identify  the 
mechanism  in  the  formation  of  triple  points  in  the  context  of  simplified  mathematical  models. 
In  (9]-[l  1]  Erpenbeck  first  used  normal  mode  analysis  on  the  linearized  Euler  system  to  study 
the  stability  of  multi-dimensional  detonation  waves.  Later,  in  [12],  Strehlow  introduced  the 
concept  of  acoustic  ray  trapping  to  study  the  formation  of  Mach  stems;  this  idea  has  been 
generalized  by  Majda  [13]  using  high  frequency  asymptotics.  With  the  rapid  advent  of 
modern  computing  capability,  another  available  avenue  in  studying  detonation  phenomenon 
is  by  direct  numerical  simulations. 

In  this  paper,  we  will  present  a  high  order  hybrid  method  in  order  to  study  tw'o  dimen- 
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sional  detonation  waves,  which  will  be  able  to  resolve  the  detailed  transverse  wave  structures 
of  multi-dimensional  detonation  waves.  The  numerical  simulations  done  in  the  paper  are  for 
an  idealized  model  of  chemical  reaction  in  which  the  reactant  species  is  in  an  irreversible, 
unimolecular  reaction  A  —*  B  with  finite  Arrhenius  reaction  rate.  It  is  evident  that  this 
model  can  not  represent  all  the  effects  of  realistic  chemical  kinetics  on  the  cellular  structures 
observed  in  lab  experiments.  However,  an  accurate  numerical  solution  of  this  simplified 
model  will  provide  a  better  understanding  of  the  physics  involved  in  the  onset  and  evolution 
of  detonation  waves  and  a  verification  of  current  mathematical  theories  on  detonation  waves. 

There  are  three  basic  characteristics  of  detonation  waves:  (1)  a  strong  precursor  detona¬ 
tion  front;  (2)  Mach  stem  configuration  of  the  “triple  points”  and  transverse  wave  structures; 
(3)  stiff  chemical  reactions.  The  flow  field  can  be  divided  into  regions  of  highly  irregular  and 
steep  gradients  near  the  detonation  front  and  regions  of  strong  but  smooth  pressure  waves. 
Near  the  shock  fronts,  strong  vorticity  fields  are  expected  from  the  roll-up  of  slip  lines.  The 
temporal  changes  of  thermodynamic  and  chemical  compositions  also  vary  dramatically  from 
region  to  region.  We  will  construct  our  numerical  schemes  according  to  these  characteris¬ 
tics  of  multi-dimensional  detonation  waves,  therefore  ,  it  is  not  surprising  that  the  resulting 
numerical  scheme  is  of  a  hybrid  type. 

The  reaction  rate  depends  on  the  flow  temperature  exponentially  through  the  Arrhenius 
relation.  Accurate  computations  of  the  flow  field  are  extremely  important  in  producing  the 
correct  chemical  reactions  and  thus  the  correct  cellular  structures.  Traditional  shock  captur¬ 
ing  schemes,  designed  to  smooth  shock  and  contact  discontinuities,  introduce  a  considerable 
amount  of  numerical  viscosity  near  those  discontinuities.  They  have  been  shown  to  have 
a  tendency  to  distort  the  real  chemical  reaction  processes.  In  [14],  the  widely  used  P.P.M. 
high-oider  Godunov  scheme  was  found  to  produce  nonphysical  weak  detonations.  Also  in 
[15],  the  ENO  finite  difference  scheme  was  shown  to  yield  wrong  detonation  speeds  in  one 
dimensional  ZND  simulations.  All  these  facts  point  out  the  importance  of  designing  numer¬ 
ical  methods  wiihoul  excessive  numerical  viscosity.  Among  the  attemts  of  simulating  two 


dimensional  detonation  waves  numerically,  Oran  et  al  has  done  a  series  of  simulations  with 
FCT  schemes  and  a  phenomenological  chemistry  model  was  used  to  solve  the  stiffness  in  the 
system,  see  for  instance  [16].  Another  early  result  was  obtained  by  Taki  and  Fujiwara  [17] 
where  a  quasi  first  order  numerical  scheme  was  used  with  a  two  front  model  suggested  in 
[18]  to  model  the  chemistry.  Most  recently,  an  improved  version  of  P.P.M.,  which  uses  the 
location  of  the  detonation  in  evaluating  the  numerical  flu.xes,  gave  improved  results  in  one- 
and  two-dimensional  detonation  simulations  [19]. 

The  work  reported  in  this  paper  is  part  of  a  larger  research  project  which  intends  to 
understand  different  aspects  of  detonation  waves,  including  the  DDT  process  and  the  effects 
of  turbulence  and  the  detonability  and  detonation  failures.  In  Section  1,  we  will  introduce  the 
governing  equations  for  reacting  flows  and  its  formulation  in  general  curvilinear  coordinates. 
In  Section  2,  we  discuss  the  hybrid  numerical  scheme  using  multi-domain  approaches,  and 
the  ENO  finite  difference  methods  and  Chebyshev  collocation  methods  and  shock  tracking 
methods  will  be  introduced.  These  numerical  techniques  will  be  used  in  the  framework 
of  multi-domain  to  fit  the  properties  of  2-D  detonation  waves.  In  Section  3,  we  consider 
the  treatment  of  interface  conditions  between  different  numerical  schemes  and  boundary 
conditions  and  the  smoothing  techniques  for  the  detonation  front.  In  Section  4,  we  validate 
the  hybrid  numerical  scheme  and  test  the  effects  of  smoothing  of  the  detonation  front  on  the 
cellular  pattern  of  detonation  waves.  Then,  we  present  the  main  results  of  this  paper,  the 
simulations  of  several  cases  of  two  dimensional  detonation  waves.  Detailed  analysis  of  the 
results  will  be  discussed.  Finally,  in  Section  5  we  give  a  conclusion  and  the  plan  for  future 
works. 
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1  Governing  Equations 


Consider  two  dimensional  detonation  waves  in  an  infinite  channel  moving  from  left  to  right 
into  unreacted  gas  mixtures  and  the  channel  is  denoted  by  Q, 


w  w 

n  =  (-00,00)X[ - y] 


(1) 


where  W  is  the  channel  width. 

The  governing  equations  for  reacting  detonation  waves  with  s  species  and  p  reaction  steps 
are  the  following  Euler  equations, 

du  ,  af(u)  ag(u) 


dT^  dx 
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dy 


=  4>(u) 
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where 


u  =  {p,pu,pv,pe,pi,---  ,psf 


u)  =  ^pu, 
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g(«) 


pu^  -h  P,  puv,  pu{€  +-),  ,(7, 


i^v.pvu. 


pv^  +  p,  pv{e  +  -),  PiU,  •  ■  • ,  P,v 


)' 


4>(u)  =  (0, 0,0,0,  a;,, 


where  (u,  v),  p,  p,  e  are  velocity  vector,  density,  pressure,  and  total  specific  internal  energy,  p, 
is  the  mass  density  of  the  ith  species.  Ui  ~  is  the  molecular  weight 

of  species  -i  and  Uij  is  the  stoichiometric  coefficient  for  the  j— species  in  the  i-th  reaction 
step,  and  ^  denotes  the  change  rate  of  the  ith  reaction  progress  variable  which  is  assumed 
to  obey  Arrhenius’  rule.  In  the  case  of  one-step  A  B  irreversible  reaction,  i.e.  s  =  p  =  1 
and 

u;  = -A>Aexp(-^)  (3) 


where  A  =  ^  is  the  mass  fraction  of  reactant.  If  we  assume  exthermoic  reaction,  the  specific 
internal  energy 
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u  ^  -f- 1>^ 
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where  Q  is  the  specific  heat  formation  and  7  =  1.2  is  the  ratio  of  specific  heats. 


All  quantities  above  have  been  non-dimensionalized  by  the  initial  states  in  the  unreacted 
gas  mixture  in  front  of  the  detonation  fronts.  They  are  given  as  follows  ”  indicates 
nondimensionalization  and  “0”  subscript  denotes  states  of  unreacted  gas  ) 
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Where  t*  =  =  xi  is  the  half  reaction  distance  which  a  mass  particle  will  travel  from 

the  detonation  front  before  half-depletion  of  the  reactant  occurs. 

Consider  a  general  curvilinear  coordinates  =  ^{x,y,T).7}  =  T){x^y,T),t  — 

t{x,y,T),  the  governing  equation  (2)  becomes 


dt^  d(  ^  dr]  ^  ’ 


io) 


where  U  =  ■y/gu, 
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lE  ^pu  +  yr^p 

% 

1 
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'I'(U)  =  {0,0,0,0, 


where  ^  —  x,,y(^  is  the  .Jacobian  of  the  mesh  transformation,  =  C  +  + 

and  +  ur/j,  -f  vrjy.  The  spatial  discretization  is  done  in  the  (C  ^/)  space  for 

equation  (5)  and  the  range  for  ^.7]  are  both  [-1,1]. 

In  most  of  the  numerical  simulations  done  in  the  paper,  for  initial  conditions  we  use 
the  Chapman-.Jouguet  ZNl)  steady  solution  which  can  l)e  obtained  by  solving  the  Hankine- 


Hugoit  equation  between  the  state  in  front  of  detonation  front  and  all  states  behind  the 
detonation  front [5]. 

2  Hybrid  Numerical  Algorithm  with  Domain  Decom¬ 
position  Technique 

In  this  section  we  describe  the  hybriri  numerical  sclu'me  for  detonation  waves.  The  conijMi- 
tational  region  is  composed  of  the  detonation  front  moving  to  the  right  and  the  rt'ar  ])iston 
boundary  and  upper  and  lower  solid  walls.  It  will  be  subdivided  using  mult  i-tiomain  tech 
nirjue  with  the  detonation  front  as  the  right  most  boundary  (see  Figure'  1).  Three  different 
numerical  techniques  will  be  applied  in  different  parts  of  the  comi)utat ional  region,  d'liey 
are 

•  Shock  tracking  algorithm  for  the  detonation  fnnit; 

•  High  Order  ENO  finite  difference  scheme  in  the  .subdomain  which  contains  t  he  re'llected 
shocks  and  contact  discontinuities  along  the  detonation  front; 

•  Chebyshev  collocation  methods  for  the  strong  vorticity  and  pressure  fields  from  the 
interaction  of  transverse  waves  along  the  detonation  front. 

We  give  a  brief  descriptions  of  each  of  the  three  numerical  techniques. 

2.1  Chebyshev  collocation  methods 

In  the  computation  domain  €  [—1,  IJXf— 1, 1],  to  approximate  any  quantity  /(^,  r;),  we 
consider  its  Chebyshev  collocation  interpolant  which  is  defined  as  (assuming  that 

polynomials  are  of  the  same  order  N  in  both  the  ^  and  rj  directions): 

0<i,i<N 

where  ^  cq  =  =  2  and  c_,  =  I,  for  1  <  j  <  A'  -  1 .  Here  ==  cos  ^ 

are  the  Chebyshev- Gauss- Labotto  points,  and  7'/v(0  —  cos(A^  ros“’{^))  is  the  (’hebyshev 
polynomial  of  the  first  kind. 
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The  derivatives  of  can  be  approximated  by  those  of  i-e. 

0<i,}<N 

=  E  E 

0<j<N  0<i<N 

where  the  inner  summation  in  the  square  bracket  can  be  evaluated  either  by  matrix-vector- 
multiplication  on  a  vectorized  machine  or  by  Fast  Fourier  Transformation.  The  latter  method 
only  involves  0{N^  log  N)  operations  for  the  computation  of  all  ?//)•  Similar  procedures 
can  be  obtained  for  fn{(k',Vi)- 

2.2  High  Order  ENO  Finite  Difference  Methods 

To  apply  the  ENO  finite  difference  method  to  (5),  the  spatial  derivative  is  discretized  by 
conservative  numerical  flux  differences: 


where  f/jj  is  the  numerical  approximation  of  (5)  using  the  method  of  lines.  In  order  to  com¬ 
pute  the  numerical  fluxes  (j- ^i,  first  the  primitive  functions  for  F,G  are  approx¬ 

imated  by  piecewise  polynomials  using  the  ENO  adaptive  stencil  idea  of  Harten  [20], and 
then  the  derivatives  of  those  polynomials  are  evaluated  at  the  edge-centered  mesh  points 


to  produce  the  numerical  fluxes.  For  technique  details  on  the  construc¬ 
tion  of  ENO  fluxes  to  the  system  of  equations  (5),  we  refer  the  reader  to  [20], [21]. 


We  need  the  eigenvalues  and  eigenvector  on  the  .Jacobian  of  fluxes  A  =  |^,B  = 
for  the  characteristic  decompositions  in  order  to  define  ENO  numerical  fluxes.  If  ^  = 
-f  denote  the  material  derivative,  then  we  have  the  following  eigenvalues  and 


eigenvectors, 
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where  a  =  is  the  sound  speed  and  Uj 
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for  A  and  = 
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primitive  variable, 


2  and  the  transform  matrix  T  =  |^  with  v  =  (p.  u,r,p.  A)  being  the 
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2.3  Tracking  Algorithm  for  Precursor  Detonation  Front 

The  detonation  front  will  l)e  represented  by  a.  continuous  curve 

,  ,  -W  W 

x  =  .r(?/,/),  <  ?/  S  I  >  b 


(11) 
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and  the  normal  of  the  front  denoted  by  n  =  (ni,ny)  points  to  the  unreacted  gas,  which  can 
be  computed  by 


Tlr 


riy  = 


v/i  + 


(13) 


-X. 


\/i  + 

By  the  Hygen’s  principal,  the  detonation  front  will  propagate  in  its  normal  direction  with 
speed 


Dr.  = 


Xt 


(14) 


Following  the  procedure  proposed  in  [23]  [22],  we  can  write  the  Rankine-Hogncit  condition 

across  the  shock  front  as  follows 

Poi^nfl  Dn)  —  Pl(^n,l 
Po{Un,C  ■  Dn)^ -j- po  =  ~  +  Pi 

^  :^^  +  iinr.,o~D,.r  +  XoQ  =  +  l{u +  X,Q 

poXo{Unfi  —  Dn)  =  PlAi(Un,l  —  Dn) 

where  “0”  denotes  the  state  in  front  of  detonation  front  and  “1”  the  states  behind  the  front. 
And  u,i  =  o  n  is  the  normal  velocity  on  the  front.  Equation  (15)  relates  the  states 

in  front  of  the  shock  and  behind  the  siiock.  To  close  the  system,  we  need  to  impose  the 
continuity  of  tangential  velocity,  i.e. 

utfi  =  ut,i.  (16) 


where  Ui  =  (— riy,  n^)  o  (ir,  v)  is  the  tangential  velocity. 

From  (15)  (16),  we  can  solve  the  quantities  (pi ,  pi,  Unj,  Ut,\,  Ai )  behind  the  shock  in  terms 
of  tho.se  in  front  of  the  shock  (/ro,  Po- *<u.o,  W(.o, -^o)-  Using  tlie  notation  in  [9],  we  have  the 
following 
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—  Ui,0 

Ai  =  Ao 


\vh(T('  w  =  -j'r,-  'iiid  K  =  ^  is  the  Marli  iitiinln’r  of  the  shock  trout  n  lati\('  to  the 

iinreartofi  ^as  and  I’o  —  —  \/l  ts  th<*  sound  spied  in  the  unreacted  mix!  me. 

In  order  to  di'rive  a  time  evi'hilion  equation  for  the  shock  front,  we  deline  tiie  tinu* 
differentiation  along  the  shock  front  ^  4  ■*■((.'/•  fttr  fixc'd  /y  €  [~\Vj'l.  11/2].  *fy 
applying  ^  on  both  sides  of  (17)  and  separating  t  he  terms  which  involves  .c, -(.*/'  f )-  obtain 
tlie  following 
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0 


54 


(7  -  1)^^  -  2 
(7  +  1)k^ 
dXo 


Nt 


2 

(7  +  1)  r' 


(Ji.t  =  Tt  ~  XytUi 


(It 


and  5  =  (1  +  {xy)^)c^,  N  =  uo  -  Xj,i’o,  T  =  Xyiio  +  i?o  and  A',  =  St 


dS 
dt  ' 


Using  ^  ~  ^  equation  (18)  can  be  rewritten  for  conservative  variable  u. 


3v  dUi 


irfu 


T^^  =  Tcx„4-Tc1 


(19) 


where  c  =  (ci,C2,C3,u,c-j,„,C4}'^,d  =  (di , dj, da.u, ^4)"^. 

Now  on  a  fixed  point  on  the  shock  front  we  consider  the  characteristic  component  of 


equation  (2)  in  the  norma!  direction  n  =  ( 


v/i  +(tv)'  '  v/l  +(xv)2 


), 


ir)u  d  z,  d  ,  . 

- L  - f  J - g  -  (J). 

dt  dn  dt^ 


(20) 


where  t  is  the  tangenlial  direction  along  the  shock  front,  f  =  Uj  f  +  »i,g.g  =  — n^f  +  n^g. 


The  eigenvalues  and  eigenvector  of  Jacobian  matrix  A  =  jj-  are 


lift  C,  (Int  “I"  C 


(21) 


hT,hT,-  A.T  (22) 

where  I,,  i  =  1,  •  •  • ,  5  are  given  in  (10)  with  n^,  =  .  1  . ,  Uy  =  . «  and  u„  =  n^u  +  iiyV. 

V>+(xy)*  V’+Uv) 

Along  the  normal  direction  of  the  shock  front,  the  (n„  4- c)-characteristic  field  approaches 
the  shock  front  from  behind,  therefore  a  compatibility  condition  ran  be  obtained  by  consid¬ 
ering  the  characteristic  combination  of  equation  (19) 

I5T-7— =  IsTcxtf -f  IsTd,  (23) 

(it 

thus  yielding  the  following  ODE  for  the  shock  speed 

X(e  =  //(u,,Uo,Xy,X(,J^(y)  (24) 
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where 


_  IsT^  -  IsTd 
hTc 

du^  _  u  J^_2_  ,  1 

dt  v/g  dt 

+  'i’(U)  is  the  residual  computed  in  (^,17)  space  behn  !  the  shock 

3  Interface  Conditions  and  Smoothing  of  Detonation 
Front 

3.1  Interfacial  Conditions  between  Subdomains  and  Boundary 
Conditions 

We  first  discuss  the  issue  of  interface  and  boundary  conditions  for  the  hybrid  numerical 
scheme.  A  correct  interface  coupling  between  the  ENO  finite  difference  method  and  spec¬ 
tral  method  is  crucial  for  the  global  stability  and  accuracy  of  the  scheme.  The  basic  idea 
behind  the  treatment  of  interface  and  boundary  condition  is  the  propagation  of  information 
along  characteristics  of  the  hyperbolic  systems  [24],  We  will  consider  the  following  situation 
separately. 

Case  1.  Interface  between  subdomains 

On  a  typical  interface  between  two  subdomains,  say  F  between  fli  and  0^  in  Figure  la, 
there  are  two  types  of  points,  i.e.  interior  point  I  and  cross  point  T. 

(a)  Interior  points  I 

Let  V  =  (P,  u,  u, /?,  A)^  be  the  primr  ive  flow  variable,  v^v’'  denote  the  solutions  com¬ 
puted  for  the  time  step  in  fl/,  fir-  Sr,i  denotes  the  normal  speed  of  the  interface  at  point 
/  with  the  normal  direction  n  =  (rtx.Uy)-  The  eigenvalues  and  eigenvectors  of  the  .lacobian 
matrix  A  are  given  in  (10)  and  (22),  and  the  corresponding  characteristic  variables  are 

?e,  =  p  —  apu,^  (2.'i) 


and 


front. 
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p-a  p 


ws  -  (IpUn  (29) 

where  denotes  an  average  state  between  and  v’',  for  instance  the  Roe-average  [28]. 

In  order  to  update  v  at  point  I  for  the  time  step  ,  we  make  the  following  correction 
on  Wi,  1  <  f  <  5  based  on  the  sign  of  the  difference  between  the  eigenvalues  and  the  normal 


speed  of  the  interface  Srj,  i.e. 


.corrected 


_  j  wl  if  A,  -  5r,/ 
~  I  if  A.  -  5r,/ 


Finally,  we  set  v*"  =  v'  =  v 


r  _  _  ^,corrected  _  cp—l  .^corrected 


Remark.  In  the  case  that  the  interface  is  the  shock  front,  'v‘  should  be  computed  by  the 
Rankine-Hogoniot  conditions  (17). 

(b)  Cross  points  between  the  wall  boundary  and  interior  interface  -  “T”  points 

In  order  to  update  the  solution  for  point  “T”  in  Figure  la  we  have  to  consider  the 
information  which  comes  from  both  subdomains  and  12^  and  also  the  role  of  the  wall. 
Characteristic  surfaces  approaching  the  point  “T”  from  several  direction  can  be  used  to 
form  a  closed  system  to  determine  this  solution  [25]  [26].  Thus,  such  an  approach  is  not 
unique  and  based  largely  on  the  experience. 

First  let  us  consider  the  characteristic  surface  which  is  tangential  to  the  interface  F  with 
normal  n  =  {nx,ny)  —  (1,0).  The  corrected  characteristic  quantities  can  be 

obtained  as  follows,  assuming  Ai,A4,A5  all  positive  (otherwise,  replace  subscript  “1”  by  “r” 
for  each  negative  values  of  A,  —  Sr,T,i  =  1,4,5), 


.corrected 


p‘  -  apul 


IC4  =  tv^ 


corrected 


—  W 


corrected 

5 


p‘  +  apix'„. 
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where  subscript  “a:”  indicates  the  normal  direction  of  the  characteristic  surface  and  =  uon. 

On  the  other  hand,  we  know  that  the  entropy  .s  remains  constant  along  the  characteristic 
direction  corresponding  to  the  eigenvalue  u„  =  uon,  we  can  correct  the  entropy  .s  as  follows 


„coTrecled 


if  Un  —  Sr,T  >  0 
otherwise 


Next,  we  consider  the  characteristic  surface  approaching  the  wall  at  point  “T”  with 
normal  n  =  {nx,ny)  =  (0,  1)  -  (top  wall)  or  (0,-1)  -  (lower  wall).  On  the  top  wall,  u„  = 
V  =  0,  so  A/  =  —  a  =  —a  and  we  have  the  first  characteristic  field  approachiiig  the  wall 

from  the  computational  domain.  Therefore,  we  can  correct  the  first  characteristic  variable 
w\  using  the  results  from  either  fli  or  i-e. 

_  ^.corrected  _  /  P'  +  OP"!,  “n  “  Sv  J  >0 


Wi  =  Wj 


p''  +  opu^  otherwise 


Finally,  we  solve  for  all  four  primitive  quantities  of  the  flow  as  follows 


u  =  0 


A  =  W4 

p  =  «  +2u;5')/4 

P  =  (-)^- 

s 


The  treatment  of  the  cross  point  on  the  lower  wall  can  be  done  similarly. 

Case  2.  Wall  Boundary  Conditions 

For  a  point  on  the  wall  which  only  belongs  to  one  subdomain  such  as  point  “B”  in  Figure 
la,  n  =  {nx:,ny)  —  (0, 1)  or  (0,-1)  for  top  and  lower  walls,  respectively.  So  u„  =  u  o  n  = 
u  =  0,  therefore  the  wall  itself  is  a  characteristic  surface  on  which  no  normal  flow  condition 
is  enforced,  i.e.  v  =  0.  Additionally,  for  the  top  wall  (lower  wall  is  treated  similarly)  we  have 
the  following  compatibility  equation  according  to  its  corresponding  outflow  characteristic 
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(41) 

(42) 

(43) 


W3  =  p~  a^p 

ai4  =  X 

W5  —  p  +  apun- 

So,  the  solution  at  point  “B”  can  be  obtained  as  follows, 

u  =  0,u  =  -  — ,p  =  105, A  -W4,p  =  (p'  -  W3)ja\  (44) 

Uy 

3.2  Smoothing  Technique  of  Detonation  Front 

In  order  to  evolve  the  detonation  front,  we  need  to  compute  the  time  derivative  of  its  normal 
speed  Dn,  which  depends  on  Xtt  in  (24).  Thus,  the  accuracy  of  the  detonation  front  depends 
on  the  residual  of  numerical  solution  at  the  detonation  front.  Numerical  experiments  show 
that  the  front  will  develop  high  frequency  numerical  instability  if  no  smoothing  is  applied  on 
the  detonation  front.  In  this  paper,  we  test  three  types  of  smoothing  on  the  detonation  front 
and  compare  the  effects  of  different  smoothing  on  the  cellular  pattern  of  detonation  waves. 
Because  detonation  waves  are  unstable  in  most  C2ise  in  the  high  frequency  range,  this  is  a 
dilemma  for  our  numerical  simulation.  On  one  hand,  we  try  to  admit  as  many  frequencies  in 
the  numerical  solution  as  possible  in  order  to  obtain  enough  nonlinear  interaction  between 
different  unstable  frequencies;  on  the  other  hand,  to  maintain  numerical  stability  some  cut¬ 
off  in  high  frequencies  is  needed  for  long  time  integration.  So,  the  best  that  we  can  hope  for 
is  to  obtain  as  fine  a  resolution  a.s  possible  in  our  numerical  simulations  while  maintaining 
the  numerical  stability. 

Smoothing  Technique  One  -  Averaging  Solution  on  the  Detonation  Front 

The  simplest  way  to  eliminate  high  frequency  on  the  shock  front  is  to  smooth  the  solution 
U]  in  H {ui,Uo^  Xy,  Xt,  Xiy)  on  the  right  hand  side  of  (24). 

Smoothing  I.  In  equation  (24)  replacing  Ui ^  by  the  simple  average  of  neighboring  solutions. 
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which  of  course  will  reduce  the  accuracy  of  the  solution  on  the  detonation  front, 

..  ,  +2u,j  +  u,,,+, 

Ui,_,  < - - - . 

It  can  be  seen  that  this  averaging  procedure  will  reduce  the  accuracy  of  the  shock  front  to 
only  first  order,  and  have  a  larger  dissipation  than  the  other  snioolhing  technicpies  suggested 
below. 


Smoothing  II  -  Derivative  Smoothing 

One  of  the  most  used  ijiethods  in  shock  capturing  scheme  to  control  oscillations  derivative 
limiting.  The  idea  was  first  introduced  in  [27]  to  construct  monotonicity  preserving  cubic 
spline. 

Let  =  l,--,nbe7i-  discrete  data  points  and  xj  <  X2  ■  •  •  <  x„.  The  cubic  spline 

n3(x)  is  a  piecewise  cubic  polynomial  which  has  a  continuous  derivative  at  nodes  x,  and 


satisfies  the  following  conditions, 


/  n3(x,)  =  Vi 

\  naixi)  =y' 


1  <  f  <  n. 


,  _  Aj-is,  + 

A._,  +  Ax. 


2  <  f  —  1 


where  ^',Ax,-  =  x^+i  —  x,.  The  derivative  at  the  end  points  is  given  by 


,  _  (2Ai  +  A2).si  —  AiSj 

“  Ax,  +  Ax2 

/  1  "I"  ^71—2)^71—1  ^Ti— l^Ti— 2 

Ax,i_j  +  Ax,i_2 


However,  the  cubic  spline  defined  as  above  is  usually  oscillatory  and  the  monotonicity 
of  the  original  data  set  will  be  lost.  It  is  suggested  in  [27]  that  the  monotonicity  can  he 
preserved  except  at  a  local  extrema  point  by  limiting  the  derivatives  ;|/'‘s. 


y, ' —  y/.m,.  = 


77tnr(7?inx(0,  if  y]  >  0 

nia.r{miii(0,y'^), if  //'  <  0. 


So  we  suggest  the  following  smoothing  proc<'dure  for  the  rietouation  front. 


SMOOTHING  II.  In  equation  (24)  change  x'y  in  H  as  follows, 

{^y)i  * —  1  <  «  <  n.  (48) 

Smoothing  Technique  III  -  High  Frequency  Spectral  Space  Smoothing 

The  third  smoothing  technique  tested  in  this  paper  applies  high  frequency  cut-off  func¬ 
tions  in  the  Fourier  transform  space  of  the  detonation  front.  Assume  that  the  discrete  data 
1  <  *  ^  ^  are  equally  distributed  and  data  (3:,,  j/j)  has  been  decomposed  into  Fourier 
modes, 

Vi  =  Yi  \  <i  <n  (49) 

where  yk  = 

We  multiply  the  Fourier  coefficient  yk  by  a  decreasing  factor  (ik  so  that  the  high  frequen¬ 
cies  will  be  decreased.  The  modified  yi  is  given  by 

^  — 1 

Vi  (t/.)^'""*"  =:  E  \<i<n.  (50) 

2 

Here  we  chose  crt  so  that  it  decays  exponentially  in  terms  of  the  wave  number, 

cr,  =  for  \k\  <  (51) 

where  the  constant  u  is  chosen  so  that  <7n  is  the  machine  zero  and  2i  is  called  the  order  of 
'  2 

the  exponential  filtering. 

4  Cellular  Structures  of  Two  Dimensional  Detonation 
Waves 

4.1  Linear  Stability  Analysis  and  2  D  Detonation  Waves 

In  [9]-[ll]  Erpenbeck  first  used  linearized  normal  mode  analysis  to  study  the  stability  of  two 
dimensional  detonation  wave.  With  complex  analysis  technique,  he  analyzed  the  unstable 
modes  of  linearized  Euler  equations  with  respect  to  the  steady  state  solution  of  plane  ZND 
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detonation  waves.  The  stability  of  the  detonation  front  llnis  tlu'  whole  tUiw  field, 

is  determined  by  the  existence  of  poles  of  its  time  Laplace  transformation  for  any 

transverse  frequency  c  =  +  It  is  found  that,  for  large  wave  numbers  (higli  frequencies), 

the  stability  of  ZND  detonation  waves  depends  on  the  quantity  (o(e)  —  i/.^(.r)  wln'ie  x  is  the 
distance  measured  away  from  the  ZND  steady  shock  front  and  co(,r)  is  the  frozen  .sound 
speed  at  location  x  and  ti{x)  is  the  flow  velocity  there.  Detonation  waves  with  an  one  ste}) 
irreversible  A  B  reaction  have  been  categorized  [1 1]  in  terms  of  their  short  wave  stabilities: 

Type  D  eo{x)  —  u^(x}  is  decreasing  as  .r  moves  away  from  the  detonation  front; 

Type  I  Co(x)  —  u'^(x)  is  increasing  as  x  moves  away  from  the  detonation  front: 

Type  M  Co(x)  — u^(x)  achieves  a  maximum  at  some  points  x*  between  x  —  0  and  x  =  — oc. 

Type  D  ha.s  been  shown  to  be  stable  for  high  frequencies,  while  Type  1  and  M  will  be 
unstable  for  bands  of  large  wavenumbers,  e,j  <  f  <  =  1,2,  •••  where  f  =  ^  and  IV’ 

is  the  channel  width.  The  latter  case  means  difficulty  in  attempting  to  simulate  detonation 
waves  through  numerical  calculations  as  it  will  never  be  possible  to  resolve  all  the  unstable 
modes  in  the  system  with  finite  a  number  of  mesh  points. 

4,2  Validation  of  Numerical  Scheme 

Computational  Mesh  Set-up  and  Initial  Conditions 

The  following  notation  will  be  used  in  all  the  computations;  ndm  -  the  number  of  sub- 
domains  in  the  subdivision  of  the  computation  domains;  W  -  channel  width;  Isk  -  number 
of  marker  points  on  the  shock  front;  (n,7n)  -  size  of  mesh  in  a  subdomain.  In  all  the  com¬ 
putations  presented  here,  we  take  ndrn  =  9  and  the  total  length  of  the  channel  to  be  loO/’* 
with  the  detonation  front  as  the  right  boundary  of  the  .solution  domain.  As  the  detonation 
propagates  and  curves,  the  interfaces  between  subdomains  will  also  travel  at  a  speed  whic  h 
is  taken  to  be  the  averaged  speed  of  tin*  curved  detonation  front  . 

In  the  first  sulxloniain.  we  use  third  order  LNO-U'  scheim's[2l)  in  order  to  resolve*  the 
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reflected  shocks  and  contact  discontinuities.  A  stretching  function  will  be  used  in  the  x- 
direction  so  that  the  mesh  will  be  clustered  toward  the  detonation  front.  The  stretcliing 
function  is  given  by 


C  =  ^{0  =  -1  +  -  sin  — 
a 


-1 


a{(  +  l) 


1  <f.r  <  1 


where  is  the  stretched  mesh  and  the  parameter  rv  determines  the  amount  of  stretching, 
we  take  a  =  0.995. 

The  right  hand  side  of  the  first  subdomain,  being  the  detonation  front,  will  be  tracked  by 
the  track  algorithm  in  Section  2.3.  Appropriate  smoothing  will  be  applied  on  the  detonation 
front  for  about  every  20  iterations;  the  exact  frequency  depends  on  the  strength  of  the  deto¬ 
nation  waves.  Chebyshev  collocation  methods  will  be  used  in  the  remaining  subdomains.  In 
order  to  ease  the  CFL  restriction  from  the  nonuniform  distribution  of  Chebyshev  collocation 
points,  another  stretching  function  (f){x,a)  is  used  to  produce  a  more  uniformly  distributed 
mesh  in  the  Chebyshev  collocation  subdomains.  We  take 


r  =  m  = 


sin  ' 
sin“*  a 


where  again  is  the  stretched  mesh  and  a  =  0.999  so  that  Chebyshev-Causs-Labotto  points 
in  ^  space  will  be  mapped  to  space  with  more  uniform  distribution. 

Finer  meshes  will  be  used  for  those  subdomains  closer  to  the  denotation  front.  A  typical 
mesh  set-up  is  given  in  Figure  lb  (only  the  first  seven  subdomains  are  shown). 

The  computation  starts  with  the  ZND  steady  state  solution  of  plane  detonation  wave. 
To  induce  the  transverse  waves,  we  introduce  perturbations  either  in  the  detonation  front  or 
in  the  flow  variables  themselves,  or  both.  The  type  of  perturbations  used  will  be  sine-like 
wave 

x(j/,0)  =:  x(y,0) -I- esin(7rj/^),  -W/2  <  y  <  W/2  (52) 


or  random  perturbations  {tranf{)). 

Reflective  solid  boundary  condition  will  be  used  in  all  computations. 
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Effect  of  Smoothing  of  the  Detonation  Front  on  the  Cvlhilnr  Patttrns 

The  procedure  of  smoothing  on  the  detonation  front  basically  introduces  t*xtra  nimierical 
viscosity  in  the  whole  scheme,  therefore  the  smaller  this  extra  viscosity  is,  the  more  reliable 
the  numerical  simulations  should  be.  This  argument  can  be  extended  to  any  numerical 
simulation  of  detonation  waves.  To  evaluate  the  effects  of  smoothing  of  the  detonation  front 
on  the  cellular  structure,  we  consider  the  ZND  detonation  wave  with  heat  formation  Q  =  50, 
activation  energy  E  =  50,  and  the  overdrive  /  =  1.6.  The  channel  width  is  taken  to  be 
IV  =  20F.  As  we  are  not. interested  in  the  detailed  structure  of  the  flow  field,  we  take  a  le,ss 
fine  mesh  -  Ef=i(n,m)  =  (50, 50)  +  (24, 50)  + (20, 50) +  (20, 50) +  (20, 40) +  (10, 20)  + (10, 10)  + 
(10,  10)  +  (10, 10).  The  number  of  marker  points  on  the  detonation  front  is  /.s/;  =  160.  The 
initial  size  of  the  subdomains  will  be  5f”  x  VK,  5f"  x  W,  5^*  x  IT,  5f”  x  V+,  lOf*  x  W,  15f'  x  IV', 
25^*  X  W,  40f’*  X  ly,  40^*  X  W.  Three  tests  are  done  to  see  the  effects  of  smoothing  on  the 
cellular  pattern.  In  all  the  numerical  results  reported  in  this  paper,  the  detonation  front  has 
traveled  at  least  20  channel  width  for  the  cellular  patterns  reported. 

Test  One  Smooth  I,  H,  and  III 

In  the  first  test  we  activate  all  three  types  of  smoothing  on  the  detonation  front.  Thus, 
strong  numerical  viscosity  will  be  produced  to  stablize  the  front.  But,  keep  in  mind,  even  in 
this  situation,  we  still  track  the  detonation  front  and  no  difference  across  the  front  is  used 
in  the  scheme.  In  Figure  2a,  we  record  the  pressure  on  the  detonation  front  for  time  i  = 
10/*  —  20/*.  A  very  regular  and  symmetric  two  cell  structure  is  produced  by  the  interaction 
of  four  different  triple  points.  This  correspond  to  a  cell  width  lOf*  -  half  the  channel  width. 

Test  Two.  Smooth  II  and  III 

In  the  second  test,  we  deactivated  Smooth  I  which  produces  th('  largest  numerical  viscosity 
among  all  three  types  of  smoothing.  In  this  ca.se,  only  one  cell  is  present  in  the  cellular  pat  tern 
(Fig  21))  which  corresponds  to  a  cell  width  20f*. 

Test  Three  Smooth  III  only 
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In  the  third  test,  we  use  only  Smoothing  III  which  uses  higli  frequency  attenuation  in 
the  Fourier  transformation  space  for  the  shock  front.  The  order  of  exponential  cut-off  in 
(51)  is  12,  thus  yielding  very  slight  numerical  viscosity  on  the  detonation  front.  The  cellular 
pattern  of  the  detonation  front  (Fig  2r)  consists  of  two  quite  irregular  cells,  a  larger  cell 
and  a  smaller  cell.  Two  major  triple  points  dominate  the  cellular  structure  along  with  tw'o 
secondary  triple  points.  A  close-up  picture  of  the  cell  pattern  is  given  in  Fig  6a  w'here  w'e 
will  further  examine  this  case  in  more  detail. 

These  tests  show  us  the  sensitivity  of  the  cellular  patterns  of  the  detonation  waves  to  the 
amount  of  numerical  viscosity  in  the  scheme.  Consequently,  we  have  to  be  very  careful  in 
applying  the  right  kind  of  algorithm  for  the  simulations  if  we  want  to  compare  the  numerical 
results  with  either  theoretical  predictions  or  experiment  results.  Presumably,  less  numerical 
viscosity  will  produce  more  reliable  cellular  patterns.  From  this  point  on,  all  numerical  results 
presented  will  use  only  the  Smooth  III  procedure  (about  every  20  iterations),  which  has  the 
smallest  amount  numerical  viscosity,  in  order  to  stabilize  the  evolution  of  the  detonation 
front.  There  have  been  several  situations  which  demonstrate  that  such  smoothing  is  necessary 
or  the  computation  will  abort  prematurely. 

Accuracy  of  Time  Integration  and  Mesh  Convergence  Studies 
1)  Comparison  of  Time  Discretizations 

The  stiffness  in  a  chemically  reacting  .system  poses  a  lot  of  difficulty  for  numerical  sim¬ 
ulations.  There  are  basically  two  issues  to  be  considered  when  choosing  a  time  integrator  , 
one  is  accuracy  and  another  is  CPU  time  efficiency.  For  the  one  reaction  system  tested  here. 
W'e  could  us  either  a  time  splitting  method  as  in  [19]  [16]  or  just  an  explicit  Runge-Kutta 
type  method.  For  the  splitting  method,  the  evolution  of  the  solution  can  be  split  into  two 
steps,  the  firsi  stej)  Iteing  an  h'uler  step  where  the  governing  equation  is  .solved  without  the 
chejnical  reaction  production  terms;  the  .secomi  stf'j)  invxjlves  only  the  reaction  term  with 
the  temperature  fx'ld  frozen  at  the  value  of  the  previous  Fuler  stej).  The  second  st('p  can  Ix' 
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explicitly  solved  here  because  only  one  reaction  is  considered.  However,  if  more  complicated 
chemistry  is  involved,  such  spl'tting  will  not  avoid  the  stiffness  problem  in  the  simulation. 

We  compare  the  results  of  the  splitting  method  (which  is  at  most  second  order)  and  the 
third  order  Runge-kutta  method  with  the  same  resolution  and  spatial  discretizations  in  the 
spatial  directions.  The  detonation  parameters  are  Q  —  50,  £  =  50,/  =  .'3;  the  mesh  sizes 
are  ELi(n,m)  =  (100,100)  +  (34,70)  +  (20,40)  +  (20,40)  +  (20,30)  +  (10,20)  +  (10. 10)  + 
(10, 10)  +  (10, 10);  and  the  channel  width  W  =  10^*.  The  initial  size  of  the  subdomains  will 
be  5t  X  W,5t  X  W,  x  x  W,  \Qt  x  W,  Ibt  x  W,  25r  x  W,  iOt  x  W,A0t  x  W.  In 

Figure  3,  we  plot  the  pressure  along  the  center  of  the  channel  at  time  T  —  ‘Mt‘.  The  solid 
line  is  the  result  obtained  by  the  third  order  Runge-Kutta  method  while  the  dots  (o)  are 
the  results  by  the  splitting  method.  We  can  see  both  results  agree  fairly  well  in  most  parts 
of  the  domain  except  that  the  splitting  method  fails  to  resolve  the  dip  in  the  pressure.  In 
the  numerical  tests  given  later,  the  third  order  Runge-Kutta  method  will  be  used  in  all  the 
computations. 

2)  Mesh  Convergence  Studies 

We  use  the  same  detonation  parameter  as  above  but  with  three  different  meshes  in  the 
spatial  direction,  the  third  order  Runge-Kutta  method  is  used  in  both  cases.  Mesh  A 
is  =  (100,100)  +  (34,70)  -I-  (20,40)  +  (20,40)  +  (20,30)  +  (10,20)  -f  (10,10)  + 

(10, 10)  +  (10, 10).  Mesh  B  is  Ef=i(n,m)  =  (120, 150)-|-(34,70)  +  (20,40)  +  (20,40)+(20,30)  + 
(10,20)-f(10, 10)  +  (10, 10)  +  (10,10).  Mesh  Bis 5:f^,(n,m)  =  (160, 200) +  (34, 70) +  (20, 40)  + 
(20, 40)  +  (20, 30)  +  (10, 20)  +  (10, 10)  +  (10, 10)  +  (10, 10).  The  initial  size  of  the  subdomains 
will  be  St  X  W,  5r  X  W,  5t  x  W,  ht  x  W,  lOr  x  W,  V2t  x  W,  25^  x  W,  AOt  x  W,  40£*  x  W. 
So  in  the  first  subdomain,  mesh  A  is  about  10  points  per  t  and  mesh  B  is  about  15  points 
per  t  and  Mesh  C  is  about  20  points  per  t.  In  Figure  4,  we  plot  the  pressure  along  the 
center  of  the  channel  ((-)  Mesh  C,  (o)  Mesh  B,  (+)  Mesh  A).  Close  agreement  can  be  seen 
among  the  results  for  all  meshes.  In  the  rest  of  the  computation,  we  will  use  at  least  the 
resolution  of  Mesh  B,  which  is  about  15  points  per  half  reaction  distance. 
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4.3  Numerical  Simulation  of  2-D  detonation  Waves 


We  present  four  cases  of  detonation  waves  using  our  high  order  liybrid  scheme,  eacli  belonging 
to  one  of  the  three  types  of  detonations  in  terms  of  the  short  wave  instability. 

CASE  I  Q  =  iO,  £  =:  50,/  =  1.2,  channel  width  W  —  IOC  -  One  cell  pal  tern 

This  is  a  case  of  Type  M  detonativ  ;i  whicli  is  unstable  for  high  frecpiencies.  The  mesh 
u.sed  in  this  case  is  =  (110,140)  +  (20.40)  +  (20,40)  +  (20,30)  +  (20,30)  + 

(10,20)  +  (10, 10)  +  (10,  10)  +  (10, 10).  The  number  of  “Marker  Points”  on  the  shock  front 
is  Isk  =  300.  The  initial  size  of  the  subdomains  will  be  5C  x  W,  5C  x  VV',  of’  x  IT,  5C  x  IT. 
IOC  X  VT,  1.5C  X  VT,  25C  x  IT,  40C  x  fT,  40C  x  W.  This  mesh  give‘s  a  resolution  of  11 
points  per  half  reaction  distance  in  the  ENO  domain.  As  a  result  of  the  high  accuracy  of 
the  Chebyshev  collocation  method,  we  find  that  very  good  resolution  of  the  flow  field  can  be 
obtained  with  far  less  points.  We  start  the  computation  with  the  ZND  steady  state  solution 
with  a  sine-like  perturbation  *  52)  on  the  shock  front  with  c  =  0.15. 

In  Figure  5a,  we  record  the  pressure  along  the  shock  front  for  time  /  =  lOt"  —  20r.  Only 
one  cell  is  produced  in  this  run  which  corresponds  to  two  triple  points  along  the  detonation 
front.  In  Figure  5b,  we  contour  five  snapshots  of  the  pressure,  temperature,  vorticity  and 
mass  fraction  at  time  T  =  30.5r,31.5C,32.5r,33.5C,34.5/*.  In  Figure  5c,  we  sketch  the 
interaction  of  the  triple  points  which  will  be  typical  for  all  the  other  later  ca.ses.  We  can 
.see  the  evolution  of  the  reflected  shock  from  the  pressure  field;  the  contact  discontinuity 
can  be  best  seen  from  the  vorticity  field.  Because  the  contour  lines  hardly  distinguish  the 
exact  position  of  contact  lines,  we  can  u.se  the  tenjperature  field  to  locate  the  position  of 
the  curving  contact  lines.  In  the  first  time  snapshot  (T  =  .30. 5t*)  of  both  Figure  5b  and  5c. 
we  see  two  reflected  shocks  (RS)  A,  B  waves  moving  toward  the  center  of  the  channel  and 
to  collide.  Following  RS  -  A,B  are  two  contact  discontinuities  respectively,  where 

the  signs  indicate  o])p()site  circulation  of  tln'se  two  contact  liiv  s  which  |)roduce  opposilt' 
vorticities.  In  snapshot  two  [T  —  31.5/*),  reflecte<i  shock  A^B  have  ('merged  from  the 
interact  ion. exchanged  dim  tions,  and  ar('  moving  away  from  each  other.  Notic«'  that  in 


the  field  of  mass  fraction,  the  layer  of  unreacted  gas  is  nuich  thick('r  Ix'hind  tlu  inci<ieul 
shock  than  the  one  behind  the  Mach  Stem.  The  chemical  reaction  in  the  layer  behuid  tlie 
incident  shock  waves  provide  the  energy  for  further  development  of  reflected  shocks  ,4.  B. 
The  contact  discontinuities  are  now  detached  from  their  triple  point  configuration.s. 

moving  downstream,  and  their  tips  are  rolled  up.  At  the  same  time  a  new  pair  of  rontacf 
lines  C2  emerge  behind  the  reflected  shock  waves  B,A  respectively.  In  snapshot  three 
{T  =  32. 5r),  RS  -  B,A  are  moving  toward  upper  and  lower  wall  respectively  and  are  ready 
to  be  reflected  away  from  the  wall  with  the  contact  lines  C2  .  C2  following  them.  In  snapshot 
three  [T  =  33. St"),  RS  B,A  have  been  reflected  away  from  the  walls  and  the  contact  lines 
C2,C~  detached  from  the  triple  point  configuration  while  an  another  new  pair  of  contact 
lines  are  created  behind  B  and  A  respectively.  In  the  last  snapshot  {T  =  34. 5t*), 

RS  -  B,A  are  ready  to  collide  again  which  finishes  one  cycle  of  the  interaction  of  these  two 
triple  points. 

Case  II  Q  =  50,  £1  =  50,/  =  1.6,  channel  width  W  =  20^“.  Two  ceils  pattern 

This  is  a  type  M  detonation  wave  which  is  unstable  for  high  frequencies,  and  we  use  the 
mesh  =  (50,250)  +  (34,70)  +  (20,40)  +  (20,30)  +  (20,30)  +  (10,20)  +  (10, 10) 

(10,10)  +  (10,10).  The  number  of  ‘Marker  Points’  on  the  shock  front  is  l.^k  =  300.  Fig¬ 
ure  6a  shows  a  two  cell  pattern  produced  by  the  trajectories  of  four  triple  points.  There 
is  a  larger  cell  with  width  approximately  10^’  (half  the  channel  width)  and  a  smaller 
one  with  width  .  Figure  6b  contains  six  snapshots  of  the  detonation  at  time  T  = 
20.25r,21.5^*,22<*,22.5<',23t*,23.5t*.  The  times  were  chosen  so  that  the  interaction  of 
the  triple  points  can  be  shown  clearly  in  the  contour  plots.  A  random  perturbation  with 
magnitude  e  =  0.3  is  used  to  perturb  the  shock  front  at  T  =  0.  Four  triple  points  are  pro¬ 
duced  (two  major  ones  and  two  secondary  ones)  and  the  cell  pattern  is  given  in  Figure  ba. 
Referring  to  the  six  arrows  and  depicted  shock  fronts  in  the  sketch  Figure  (ic  which  corre¬ 
sponds  roughly  the  six  time  snapshots  in  Figure  61)  .  In  the  first  time  snapsiiot  T  —  20.25/', 
in  the  lower  middle  part  of  the  channel,  two  triple  points  (\  D  an’  moving  away  from  each 


other  just  after  collision  and  .4  has  been  reflected  from  the  iijjjxT  wall  and  is  ahuul  tu  collide 
with  triple  point  B.  In  snapshot  two  (7’  =  21. of),  in  tin*  upper  j)art  of  the  channel,  two 
smaller  triple  points  .4,  B  collide  aiul  separate  and,  in  the  middle  ])art  of  the  i  haimel  triple 
point  C  is  moving  upward  to  collide  with  triple  point  A.  Near  the  lower  wall,  triph-  point 
D  is  reflecting  away  from  the  lower  wall.  In  snapshot  three  {T  —  22t‘),  triple  point  B  is 
approaching  the  upper  wall,  triple  point.s  ,4  and  ('  collide,  and  triide  ])oint  D  keeps  moving 
up  away  from  the  lower  wall.  In  snapshot  four  [T  —  22. 5r),  triple  point  B  reflects  away 
from  the  upper  wall  and  is  about  to  collide  with  triple  point  C,  while  triple  points  .4  and 
D  are  about  to  collide  with  each  other.  In  snapshot  five  {T  =  23/*),  in  the  upper  part  of 
the  channel,  triple  points  B  and  C  collide  and  triple  points  A  and  D  approach  each  other. 
Finally,  in  snapshot  six  (T  =  23.5/’),  triple  points  B  and  C  finish  the  collision  and  exchange 
directions  while  triple  points  A  and  D  collide. 

Also  notice  that  in  Snapshots  4  and  5,  the  two  pressure  waves  from  the  reflected  shocks 
intersect  with  each  other  before  the  collision  of  the  triple  points  happcuis  along  the  detonation 
front  (Snapshot  6).  Such  interaction  will  cause  sudden  reaction  of  any  unreacted  gas  in  the 
interior  region  and  produce  so-called  “explosions  in  explosions’. 

Case  III  Q  —  50,  E  =  50,/  =  3,  channel  width  W  —  107*.  One  cell  pattern  and  chaotic 
flow  fields. 

This  case  repre.sents  strong  heat  release,  large  overdrive,  and  belongs  to  type  I  wdiicli  is 
again  unstable  for  a  range  of  high  frequencies.  We  use  a  mesh  J2^--i{n,7n)  =  (120,  150)  -f 
(34,  70)  +  (20, 40) -f  (20, 30) -+-(20, 30) -|-(  10, 20) +  (10,  !0)  +  (10, 10)  +  (10, 10);  and  the  number 
of  ‘.Marker  points’  on  the  shock  front  is  Isk  —  300.  The  initial  size  of  the  ,sul>domain.s  will  be 
57*  X  W,  ryt  X  W,  87*  X  \\\  5r  x  IT,  !  07*  x  W,  1 27*  x  W,  257*  x  W,  407*  x  H',  407*  x  IV,  Only  t  wo 
triph'  |)oints  ar('  produced  in  this  case;  an  oiU'  cc'll  pattern  of  the  detonation  front  is  given 
in  Figure  7(a).  In  P'igure  7(1)),  we  have  six  time  snapshots  of  the  pressur*'.  tem|)erature. 
vorticity  and  mass  fraction.  Large  vorticities  are  eenerat('<l  behiini  the  detonation  front  and 
the  flow  field  becomes  quite  chaotic.  4  here  are  oidy  two  triple  points  along  tin'  detonation 


front  which  produce  a  one  cell  pattern  for  the  cellular  structure  with  cell  size  W  =  lOf  (see 
Figure  7a). 

5  Conclusion 

We  have  developed  a  high  order  numerical  scheme  which  is  suitable  for  computing  detailed 
transverse  wave  structures  of  two  dimensional  detonation  waves.  The  numerical  algorithm 
uses  a  multi-domain  approach  so  that  different  numerical  techniques  can  be  applied  for 
different  components  of  detonation  waves.  The  propagation  of  waves  across  the  interfaces  of 
subdomains  have  been  very  smooth  and  the  order  of  accuracy  of  the  whole  numerical  scheme 
is  only  limited  by  the  accuracy  of  the  time  integrator.  Tracking  of  the  detonation  front  avoids 
differences  across  the  detonation  front,  thus  avoiding  excessive  numerical  viscosity  in  shock 
capturing  schemes.  The  high  resolution  of  the  Chebyshev  collocation  method  enables  us  to 
use  far  less  grid  points  in  most  of  the  solution  domain  and  yields  great  savings  in  the  total 
CPU  cost.  The  potential  for  using  higher  ENO  finite  difference  schemes  in  the  subdomain 
which  contains  the  reflected  shocks  and  contact  discontinuities  can  be  further  exploited. 

We  have  shown  that  the  cellular  pattern  of  the  detonation  waves  is  affected  by  the 
accuracy  of  the  detonation  front  and  the  amount  of  numerical  viscosity,  especially  the  amount 
of  viscosity  involved  in  the  time  evolution  of  the  detonation  front  by  the  numerical  scheme. 
We  believe  that  this  point  should  be  well  taken  in  the  further  investigation  of  detonation 
waves  in  order  to  have  meaningful  comparisons  with  experiment  results. 

We  have  studied  several  cases  of  detonation  waves  with  specific  ratio  7  =  1.2,  from  small 
heat  release  (Case  I)  to  large  release  (Case  II  and  III)  and  small  overdrives  (Case  I)  to 
large  overdrives  (Case  II,  III).  The  numerical  results  successfully  reproduced  the  onset  and 
evolution  of  the  transverse  wave  structures.  The  contact  lines  within  triple  points  create 
large  vorticity  fields  behind  the  detonation  front  which  will  distort  and  interact  with  the 
detonation  front.  The  contact  discontinuities  from  the  triple  points  after  their  collisions 
convect  downstream  and  generate  vorticity  downstream.  Further  work  will  be  done  by  using 
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more  realistic  chemistry  models  so  that  comparison  with  experinx-nt  results  will  lx*  pt^ssible. 
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Figure  1  (a)  Multidomain  set-up  for  the  hybrid  numerical  scheme  for  detonation  waves;  (b) 
Mesh  structures  for  2-D  detonation  waves 
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Figure  2  Cellular  pattern  for  Q  =  50,  E  =  50,/  = 
detonation  front  (a)  Smoothing  I,  II  and  III,  (b)  Smc 
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Figure  3  Pressures  along  the  center  of  channel  obtained  by  time  splitting  (0)  and  third 
order  Runge-kutta  (-). 


Figure  4  Mesh  convergence  studies:  Pressure  along  the  renter  of  channel  with  three  meshes 
(-)  20  points  per  t,  (o)  15  points  per  t,  (+)  10  points  per  t. 
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Figure  5  (Jase  I,  Q  =  10,  £’  =  50,/  =  1.2,  VF  =  10  (a)  one  cell  pattern,  (b)  five  snapshots 
of  pressure,  temperature,  vortirity,  mass  frartion  of  reartant  (from  top  to  bottom)  at  T  - 
•■10.5<*,  31.5f*,  32.5f*,  34.r)r,  (c)  sketeli  of  the  interac  tion  of  the  tri]>le  |)oiuts. 
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Figure  6  Case  II,  Q  =  50,  E  =  50,/  =  1.6,  VF  =  20  (a)  two  cell  pattern,  (b)  six  snapshots 
of  pressure,  temperature,  vorticity,  mass  fraction  of  reactant  (from  top  to  bottom)  at  T  = 
20.25r,21.5r,22r,22.5r,23r,23..5r.  (c)  Tracks  of  the  triple  points. 
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Figure  6b 
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Figure  7  Ctise  III,  Q  =  50,  E  =  50,/  =  Z,W  —  10  (a)  one  cell  pattern,  (b)  six  snapshots 
of  pressure,  temperature,  vorticity,  meiss  fraction  of  reactant  (from  top  to  bottom)  at  T  = 
33r ,  33.5r ,  34r ,  34.5r,  35r ,  35.5r . 
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